Behavioral analysis has become one of the most important advancements in the movement ecology literature. Due to the importance of the internal state in individual movement decisions, analyses such as these offer an opportunity to understand the motivations underlying space-use more clearly than most of the broader scale analyses we’ve seen thus far. We are going to begin by a simple example of First Passage Time (FPT), then consider more complex methods such as Hidden Markov Models (HMM) and Behavioral Change Point Analysis (BCPA).
library(tidyverse)
library(adehabitatLT)
library(sf)
library(moveHMM)
library(bcpa)
#load albatross data
albatross <- read_tsv("data_files/Study2911040", comment="##", col_types = "Tnncc") %>%
as.data.frame() %>%
na.omit()
#add UTM coordinates
coords <- sp::SpatialPoints(data.frame(albatross$location_long, albatross$location_lat), proj4string=CRS("+proj=longlat +ellps=WGS84"))
UTM_coords <- spTransform(coords, CRS("+proj=utm +south +zone=16 +ellps=WGS84")) # reproject to UTM (meters)
albatross$utm_long <- UTM_coords@coords[,1]
albatross$utm_lat <- UTM_coords@coords[,2]
First passage time (FPT) measures the length of time an animal spends within an circle of a given radius and is commonly used as a proxy for area-restricted search behavior (Fauchald & Tveraa 2003). We first calculate the FPTs for a range of radii, and plot the variance in FPT vs. radius to identify a ‘characteristic scale’ of restricted search behaviour where the variance is maximized. FPT is easily calculate in the adehabitatLT package, so first we need to convert our data into an ltraj object. Let’s start with an FPT analysis on our Galapagos Albratoss tracking data.
# create 'ltraj' object
alba_ltraj <- as.ltraj(xy=cbind(albatross$utm_long, albatross$utm_lat),
date=albatross$timestamp,
id=as.factor(albatross$individual_id),
proj4string = CRS("+proj=utm +south +zone=16 +ellps=WGS84"))
# Calculate fpt in hours for one individual (#2911074) in radii of 0-400 km in 10-km increments
albatross_fpt <- fpt(alba_ltraj[10], radii = seq(0, 400000, 10000), units="hours")
# Calculate the variance for log(FPT) for radii and find the peak radius
fpt_var<-varlogfpt(albatross_fpt)
There appears to be a peak at very small spatial scales (1 km) and a secondary peak at large spatial scales (about 120 km). Let’s investigate the behavioral patterns at large spatial scales.
# plot fpt over time
plot(albatross_fpt, scale=120000, xlab=c("date")) #scale = radius to plot at
# add fpt at characteristic scale to tracking data
albatross10 <- albatross %>% filter(individual_id=="2911074")
albatross10$fpt <- albatross_fpt[[1]]$r12
# plot track by fpt
ggplot(albatross10) +
geom_point(aes(x=location_long, y=location_lat, color=fpt)) +
scale_color_gradientn(colors=rev(topo.colors(10))) +
theme_bw()
Here we can pick out based on FPT high-use areas around the albatross’s breeding area on the Galapagos Islands to the East, and foraging areas during the non-breeding season off of mainland South America to the East.
The next method we’re going to take a look at is the behavioral change point analysis (BCPA), which looks for the points in a time series during which there are notable shifts. In our case, we will be applying the method to a movement trajectory to see where an animal may transition between behavioral states, but technically change point analyses can be performed on any time series data (e.g., fluctuating stock values over time or carbon dioxide concentration in the atmosphere over time).
Just as with all other packages, bcpa has its own data format that it prefers, so we will use the bcpa::MakeTrack command to translate a set of 100 coordinates (from our 1326 point albatross path, for the sake of readability in the outputs) into a usable format:
alba_bpca <- bcpa::MakeTrack(albatross10$location_long,albatross10$location_lat,albatross10$timestamp)
plot(alba_bpca)
To obtain the step length and turning angles, use the bcpa::GetVT command, which decomposes the data into single steps and calculates all the statistics:
albatross.VT <- GetVT(alba_bpca)
head(albatross.VT)
## Z.start Z.end S Phi Theta
## 2 -89.62097-1.38952i -89.62086-1.38949i 1.055716e-04 0.2254499 2.967982
## 3 -89.62086-1.38949i -89.62098-1.38952i 1.225402e-04 -2.9336341 -3.159084
## 4 -89.62098-1.38952i -89.62092-1.38953i 6.858608e-05 -0.1862436 2.747391
## 5 -89.62092-1.38953i -89.62101-1.38951i 9.487602e-05 2.8585153 3.044759
## 6 -89.62101-1.38951i -89.62095-1.38952i 5.426389e-05 -0.2628578 -3.121373
## 7 -89.62095-1.38952i -89.62098-1.38952i 2.901672e-05 2.9719114 3.234769
## T.start T.end T.mid dT V T.POSIX
## 2 0.05416667 1.532222 0.7931942 1.478055 7.142605e-05 2008-06-23 18:46:19
## 3 1.53222167 3.034999 2.2836106 1.502778 8.154246e-05 2008-06-23 20:15:44
## 4 3.03499944 4.527778 3.7813886 1.492778 4.594525e-05 2008-06-23 21:45:36
## 5 4.52777778 6.036944 5.2823608 1.509166 6.286652e-05 2008-06-23 23:15:40
## 6 6.03694389 7.533333 6.7851386 1.496389 3.626322e-05 2008-06-24 00:45:50
## 7 7.53333333 9.035833 8.2845832 1.502500 1.931230e-05 2008-06-24 02:15:48
The essence of a change point analysis is a sweep across a time series in search of breaks. This sweep can be conducted in a number of ways, but we will focus here on the window sweep, whereby we identify an appropriate windowsize and sensitivity (K) and then the algorithm searches across the time series in search of break points. One can also input a function as the second argument (it can represent any combination of the elements of our albatross.VT dataframe), to serve as a response variable. In this case, we will define a very simple function persistence of movement in a given direction (i.e. ‘persistence velocity’) because we dont really have any a priori conception of what exactly causes change points in this path.
albatross.ws <- WindowSweep(albatross.VT, "V*cos(Theta)", windowsize=100, progress=FALSE, K=2)
The object that is returned by this function (which takes a little while to run, hence our reduction of the dataset to a smaller length) is a ws data frame whose final column indicates proposed break points should be and the parameter values associated with before and after those break point.
head(albatross.ws$ws)
## Model LL bic mu1 s1 rho1 mu2
## 1 7 720.7501 -1409.194 -2.458589e-05 3.673670e-05 1.0837768 0.11439631
## 2 7 720.7322 -1409.159 -2.449844e-05 3.663456e-05 0.9758211 0.11722744
## 3 7 712.8886 -1393.471 -2.377652e-05 3.629145e-05 0.9474522 0.11278957
## 4 7 705.5419 -1378.778 -2.353749e-05 3.646370e-05 0.9346491 0.10748713
## 5 7 697.9346 -1363.563 -2.303054e-05 3.642509e-05 0.9205235 0.10344379
## 6 7 690.3318 -1348.358 -2.285653e-05 3.663489e-05 0.9218635 0.09985845
## s2 rho2 Break.bb.time
## 1 0.1337935 7.588469 119.2833
## 2 0.1318655 7.255237 120.7824
## 3 0.1303602 7.361823 120.7824
## 4 0.1298768 7.631846 120.7824
## 5 0.1285572 7.821147 120.7824
## 6 0.1271208 7.999981 120.7824
We can take a look at these suggested breakpoints by looking at the smoothed plot (i.e., the summary in which all the windows are averaged to obtain the “smooth” model). In this plot, the vertical lines represent the significant change points, the width of the lines is proportional to the number of time that change point was selected.
plot(albatross.ws, type="smooth", mu.where="topright")
That doesn’t offer the clearest picture. We can see that there are many change points that have some support. We could, however, add a threshold parameter, which indicates how many of the windows that were swept over the data must have selected a particular changepoint for it to be considered significant. Here, we will use 20 and see what it looks like:
plot(albatross.ws, type="smooth", threshold=20, mu.where="topright")
This reduces our number of change points down to a more reasonable 8, and all of them appear to signify reasonable shifts in our response variable (which combines velocity and angle).
An alternative way to search for change points is to use the ‘flat’ rather than ‘smooth’ method. This analysis first selects changepoints that it deems significant by clustering neighboring change points, and then estimates a homogeneous behavior that occurs between those changepoints.
plot(albatross.ws, type="flat", mu.where="topright")
Once again, if we don’t set an equivalent to the threshold parameter (in the case of the ‘flat’ approach, its called clusterwidth), we get quite a few change points. If we set this parameter to 20, we get the following:
plot(albatross.ws, type="flat", clusterwidth=20, mu.where="topright")
This fairly conservative approach results in 19 significant change points in our time series. A summary of these change points can be obtained using the bcpa::ChangePointSummary command. Here ‘mu.hat’, ‘s.hat’, and ‘rho.hat’ are the estimated mean, standard deviation, and autocorrelation of our response variable (velocity persistence) in each phase.
summary <- ChangePointSummary(albatross.ws, clusterwidth=20)
head(summary$phases)
## t.cut mu.hat s.hat rho.hat t0 t1
## 1 (-0.207,120] -2.458589e-05 3.673670e-05 1.08377685 -0.2068058 120.0707
## 2 (120,237] 9.173512e-02 1.286393e-01 4.86922976 120.0706708 236.8816
## 3 (237,392] 5.350859e-02 9.301482e-02 1.89849136 236.8815555 392.3500
## 4 (392,496] 8.122618e-03 4.498000e-02 0.03731658 392.3500292 495.9869
## 5 (496,560] 1.855522e-01 1.919126e-01 2.85542760 495.9868900 559.6516
## 6 (560,692] -1.945615e-05 7.206794e-05 1.08645573 559.6516260 692.4168
## interval
## 1 120.27748
## 2 116.81088
## 3 155.46847
## 4 103.63686
## 5 63.66474
## 6 132.76521
This summmary suggests eight phases. We can also visualize the path itself with the associated change points using the bcpa::PathPlot command or the bcpa::PhasePlot command:
par(mfrow=c(1,2))
PathPlot(alba_bpca, albatross.ws, type="flat", clusterwidth = 20, main="Flat BCPA", plotlegend=TRUE)
PathPlot(alba_bpca, albatross.ws, type="smooth", main="Smooth BCPA", plotlegend=TRUE)
We see similar patterns as in our previous methods, with periods of high direction persistence during transiting between breeding/foraging areas.